Landslide detection and inventory updating using the time-series InSAR approach along the Karakoram Highway, Northern Pakistan

Karakoram Highway (KKH) is frequently disrupted by geological hazards mainly landslides which pose a serious threat to its normal operation. Using documented inventory, optical imagery interpretation, and frequency-area statistics, the features of slope failure, the spatial distribution, and their link to numerous contributing factors have all been effectively explored along the KKH. An updated inventory for the area was recreated using the interferometric synthetic aperture radar (InSAR) persistent scatterer (PS) technology to further investigate millimetre-accurate measurements of slope deformation (Vslope). Utilizing the PS approach, Sentinel-1 data from Jan 2018 to Jan 2022 were processed by which we obtained a deformation rate (VSlope) that varies between 0 and 364 mm/year. A total number of 234 landslides were cited from the literature and classified while 29 new potential landslides were detected and several pre-existing landslides were redefined by the InSAR approach, which was incorporated to generate an updated landslide susceptibility model with 86.6% of prediction precision in the area under curve method. As previous studies done by applying the InSAR technique incorporated a short span temporally and they missed some highly deforming zones like Budalas and Khanabad landslides, contain mean velocities > 50 mm/yr, which we studied individually in this work. In this study, a comprehensive application of the InSAR technique to assessing its performance in detecting and analysing landslides has been applied. The deformation velocity (Vslope) model shows high displacement in some regions, which needed to be further investigated by geoscientists, and the updated developed landslide inventory and susceptibility map can be used for land use planning and landslide mitigation strategies.

www.nature.com/scientificreports/ fluvial landforms. The MKT and Karakorum Fault (KF) cause brittle deformation, which are the major tectonic features in this region, due to which rock masses are severely fractured and jointed 28 . Lithology plays an important role in triggering landslides. The area is characterized by fractured and weathered rock masses possessing diverse igneous, metamorphic, and sedimentary rocks (Fig. 2). The Baltit Group, Chalt schists, Quaternary sediments, Gujhal dolomite, kilk formation, and deformed Misgar slates are the most prominent local lithologies, and they are all tectonically affected and responsible for slope destabilization alongside the highway 7 . The Southern Karakoram metamorphic Complex (SKm), the Hunza Plutonic Unit (HPU), the Shaksgam Formation (SF), and Quaternary (Q) Deposits comprise the region's geology 26 . Paragneises with interbedded pelites and amphibolite constitute the SKm. Permian massive limestones are part of the Shaksgam Formation, a region of the northern Karakoram landscape. Plagioclase, quartz, biotite, and hornblende are found in the HPU, a portion of the Karakoram batholith. Lithologies of different ages are exposed along the KKH have been weathered and weaken by seismic, hydro-climatological, and anthropogenic activities leads to huge landslides and land deformation in the area.

Methods and materials
Datasets. In this study, 4 years (Jan. 2018 to Jan. 2022) imagery of the C-band Sentinel-1 SAR dataset was downloaded from Alaska Satellite Facility (ASF) online system (https:// search. asf. alaska. edu), which included scenes in ascending and descending paths as shown in Table 1. For the processing of SAR data for the ground deformation investigation and the InSAR time-series analysis, MATLAB and SARPROZ platforms were employed. The Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) of the 30-m resolution was taken from the USGS database (https:// earth explo rer. usgs. gov) to extract the topographic parameters in ArcGIS environment. To create landslide inventory Google Earth imagery was incorporated with the existing data which was validated during fieldwork. Sentinel-2 optical sensor data of 10-m resolution was also received from the USGS Earth Explorer database system for supervised landcover mapping in R-studio package. Rainfall and landslide events are directly proportion-related, so to evaluate it annual precipitation data was accessed via the Climate Hazards Group Infrared Precipitation with Station data (CHIRPS) online database system (https:// www. chc. ucsb. edu/ data/ chirps). Geology maps and fault lines were extracted by the modification of maps designed by 12,26 in ArcMAP 10.5 software. The methodological flowchart has shown in Fig. 3.

Landslide inventory.
An inventory map has information about the location, date of occurrence, and types of landslides that have left discernible traces in an area 29 . It is a record of all kinds of past landslide events in the area. Landslide inventory maps are developed for a variety of reasons, such as identifying the location and type of landslides in a given area, illustrating the effects of a single landslide-triggering event, such as an earthquake, an intense rainfall event, or a rapid snowmelt event, highlighting the abundance of mass movements, calculating the frequency-area statistics of slope failures, and supplying pertinent data to build landslide susceptibility models or hazard models 30 . In this study, 234 landslides were mapped using past studies 7,8,12,28,31,32 , optical   (Fig. 4). This inventory map was applied to validate the InSAR-based detection of landslides in the area. After InSAR processing, some new landslides were mapped according to high deformation velocity (V slope ) and the inventory map was updated accordingly.
Landside conditioning parameters. The terrain, geology, tectonic features, weather, land cover, and anthropogenic variables in the area greatly influence the spatial distribution and intensity of the landslides. Since there is no universal agreement on which triggering parameters should be included in modelling for landslides since they are very difficult-to-understand natural phenomena, it is important to include controlling factors that are pertinent to the study area and have access to accurate data 24 . According to previous studies slope, relief, fault lines, geological setting, precipitation, and barrenness in the area are the most responsible factors in landslide triggering in the study area. In this study, these six factors (Fig. 6) were taken under investigation to check the correlation with the inventory data. The spatial relationship between a landslide occurrence location and each landslide-conditioning factor was derived using the frequency ratio (FR) model. It is crucial to construct a comprehensive landslide inventory consistently over a sufficiently broad spatial extent and a long period of time. Moreover, by providing spatial-temporal coverage, several potential landslide conditioning factors can be statistically associated 33 . The frequency ratio measures the proportion of the study region where landslides happened as well as the probability of a landslide occurring versus not occurring for a particular attribute 34 . FR was calculated using the following equation; where dA is the area of landslides in the given class, dB is the area of the class; ∑dA is the total sum of the area of the landslide in the entire study area, ∑dB is the sum of the entire study area. If the value of FR is greater than 1, it shows a high correlation, and less than 1 means a very low correlation with that class.
PSInSAR processing. The spaceborne interferometric radar technique is a potent instrument for spotting ground motions on the Earth. Due to the excellent cost-benefit ratio, non-invasiveness, large area coverage, and  www.nature.com/scientificreports/ high precision of satellite data analysis, mapping geomorphologic processes and monitoring slope instability can both considerably benefit from it. to detect landslides PSInSAR technique is effective and applied by many researchers in different parts of the world [20][21][22][23]35,36 . The Permanent Scatterers (PS) approach uses co-registered, multi-temporal synthetic aperture radar (SAR) imagery (at least 15) to analyse a backscattered signal to find out highly reflecting ground features that are stable from an electromagnetic aspect 20 . The spatiotemporal baseline of PSInSAR have been calculated for descending path and ascending path is shown in Fig. 7. To create twodimensional complex image maps of the surface with dimensions ranging from slant range to line of sight (LOS) range, SAR systems, which are active, use microwave lights and save the electromagnetic echoes reflected from the surface. The displacement time series and displacement rates of each stable point (PS) can be calculated along the SAR LOS concerning a reference point that is supposed to be stable when a significant number of independent radial light and radial phase stable points (PS) exist within a radar scene and enough radar acquisitions have been collected 37 . The displacements recorded for each PS are calculated using a stable ground point with known coordinates as the reference point. Different sensors having the same wavelength collect data throughout time and can be used in a multi-temporal analysis of ground deformation. This indicates that each PS dataset's ground deformation velocity is measured along a separate LOS 38 . PSI focuses on stable-point scatterers (having the same wavelength) that are unaffected by speckles and provide a better signal for analyzing data. Permanent Scatterers (PS) are stable points that return the signal with different durations of time when looking in the same location, which allow for long-term ground subsidence calculation. The interferometric phase (Ø Int ) of a SAR signal of wavelength λ between two different images can be expressed as: In the above Eq. (2), Ø topography is represents phase change due to the relief, Ø Movement is the terrain motion in the difference between the two different images taken at a different time. Whereas Ø Noise reveals the phase noise which includes other noise components as well and the last is the Ø Atmosphere which represents the phase component due to atmospheric disturbances.
In this research work, Sentinel-1 sensor Single Look Complex (SLC) having polarization VV data was used which consists of a constellation of two polar-orbiting satellites that operate day and night and use C-band (2) ∅Int = ∅topography + ∅Movement + ∅Noise + ∅Atmosphere where D A represents the amplitude dispersion, m A is the mean deviation of amplitude in time, and σ A is the standard deviation of amplitude in time.
These strict criteria are only satisfied by a small subset of points, yet it is essential for accurate APS calculations. This mountainous area utilized only those point targets, which had amplitude dispersion values less than 0.25, the selection of pixels below the 0.25 threshold of ASI allowed selection points of high ASI. This also ensures the selection of only those PS points owning minimum decorrelation noise.
A reference network must be created by linking the PSs using Delaunay triangulation after choosing the first order PS. Each edge's differential residual topographic error (RTE) and differential deformation velocity are calculated. The APS is then computed from the phase residuals by graph inversion after the estimated linear model www.nature.com/scientificreports/ (linear displacement velocities and residual height) is eliminated. In Second order PS selection, the requirement was lower (ASI > 0.6) to produce a PS collection that was denser, where 213,601 points for descending and 298,827 points for the ascending track were selected. Then, using the same parameters and the same reference point as for APS estimation, the final procedure with APS removal was carried out. The V LOS was calculated but it does not represent actual ground target displacements. By integrating the geometric relationship between the SAR sensor and the terrain, it is possible to utilize the directions of velocities along the steepest slope gradient (V slope ) to get over V LOS ' limits. V slope is thought to represent the direction in which actual deformations caused by potential slope failure will occur most frequently 22 . The rate of deformation in the LOS direction is insufficient for reflecting the true deformation of the slope in mountainous areas 3 . The conversion for V LOS to V slope was done using the below Eq. 39 . In the final step, to assess potential landslides deformation velocity was classified. The criteria for threshold selection have been done by different researchers according to the nature of the study area and the purpose of the study by incorporating different InSAR techniques. The degree of activity of landslides can be assessed in order to calculate the average V slope of landslides that contain a significant number of coherence thresholds (CTs) and define the V slope stability threshold statistically 40 . The matrix technique uses multi-temporal InSAR datasets as indicators of the activity and intensity of landslide processes 41 . T0 identify landslide occurrence, a threshold of − 20 mm/y along the slope direction was established 31 and some other methods have been applied in different studies. Depending on the lithological properties, failure mechanisms, sensor measurement accuracy, and the aims of the inquiry, further interpretation of the observed pattern can lead to the identification of the displacement rate threshold to identify potential landslides and the statistical factors such as the standard deviation of displacement rates are frequently used to inform these criteria 42 . In our study we were aimed to map the high potential landslides, therefore > 25 mm/yr of V slope threshold were implied.
Landslide susceptibility mapping. The mapping of locations with an equal likelihood of experiencing landslides during a given timeframe is known as landslide susceptibility. The assessment of the terrain's susceptibility for a slope failure, in which the susceptibility of the terrain for a hazardous process expresses the likelihood that such a phenomenon occurs under the specified terrain conditions or parameters, and the estimation of the likelihood of a triggering event constitute a landslide hazard zonation 43 . Analysis of the susceptibility from landslides can help design policies for land use planning and offer helpful information for preventing catastrophic damage. By identifying the factors that affect landslides, estimating the relative contribution of slope failure-causing factors, establishing a relationship between the factors and landslides, and making predictions about future landslide hazards based on that relationship, the analysis is used to understand the factors that affect landslides. For landslide susceptibility mapping, inventory data is most significantly related to mapping of all landslides and spatial accuracy, so InSAR-based updated inventory was used to develop a refined susceptibility map for the area with the incorporation of the FR method in R-studio (Eq. 1).

Results
Deformation detection. The Earth's surface deformation from Jan 2018 to Jan 2022 was calculated by using 0.7 as the temporal threshold for coherence. The displacement velocity along the LOS (V LOS ) is between − 346 and 376 mm/yr for descending track and − 364-364 mm/yr for ascending mode having a total number of 512,428 PS points in the area (Table 1 and Fig. 8). A large number of PS points were found in barren land as compared to vegetated areas.
Landslide mapped based on deformation velocity (V slope ). In this study, the Budalas area was spotted as highly deformed where the mean deformation velocity (Vslope) is > 25 mm/yr (Fig. 9). The displacement movement is towards the south and the slope gradient varies from 15 to 60 degrees. It is a complex landslide, that contains steady velocity from top to bottom (Fig. 9) but the upper part exhibits the highest deformation. The lithology belongs to the Southern Karakoram Metamorphic complex (SKm) which has the most number of landslides in the area 8 , and the area is located near MKT and other small active faults. www.nature.com/scientificreports/ Khanabad landslide has complex nature and is the most deforming zone in this study possessing a mean Vslope > 50 mm/yr, which was also detected by 12 . The displacement movement is southward and has high deformation velocity at the top. The gradient steepness is from 20 to 65 degrees and is in SKm formation. The lower part has a rockfall slide as seen in optical images (Fig. 10) and the upper side has scree moments.
Mayoon area has a complex landslide having an average deformation displacement (V slope ) > 20 mm/yr (Fig. 11). This landslide has several parallel cracks (0.1-5 m opening) 28 and slides towards the south. The area has 15 to 45 degrees of slope gradient comes in the locality of the Yasin group. In 1976 first landslide was triggered in the area ruining agricultural land to a small extent. The infrastructure was damaged in 2011 when a second slope failure occurred on the eastern end of the escarpment. Twenty families had to be evacuated from the area in 2012 because of a second, comparatively smaller-scale catastrophe 28 .

Correlation between conditioning parameters and updated inventory. The landslide inventory
map for the area was converted into a raster form to calculate the number of pixels in different classes of conditioning parameters ( Table 2). The assessment showed that, more than 80% of landslides in < 45 degrees of slope (Fig. 12, Table 2). The gradient classes from 15 to 60 degrees are positively correlated with landslide activities because the terraces of the Hunza river accumulated with ancient landslides, which are unstable in the area. According to the statistical estimation the elevation areas from 1500 to 3600 m have around 90% of landslides ( Table 2). The class 1800-2600 m and 2600-3600 m have the FR value of 2.56 and 2.05 for landslide, which shows a high correlation. The sloped body having less than 3600 m of relief has seen extensive weathering, and gravity, in addition to fluvial and runoff erosion, is the primary driver of slope collapse. The zones that are higher than 3600 m of elevation are found in exposed bedrock locations that are less impacted by river and rainfall erosion, and hence low numbers of landslides occur.
The area has active tectonic nature and experienced many earthquakes in the past, so on the premises of fault lines several landslides were mapped. Based on the presumption that slopes near a fault are more likely to experience deformation, the distances from various faults were quantitatively estimated, which shows there are many densely populated landslides near active faults (Fig. 12). The strength and deformability of the bedrock are directly impacted by the presence of faults. A distance of less than 5 km is positively correlated with landslides in the area (  Fig. 12). In short, deformation hazards mostly occur in sedimentary and metamorphic rock along the KKH.
KKH is a semiarid zone have less annual precipitation (Table 2). Statistical analysis shows those areas which experience high rainfall possess a large number of slope failure events ( Table 2). There are around 95% of   13,14) were mapped in non-vegetated zones. PSInSAR method has limitations in vegetation-covered areas, but more than 60 percent of the land is covered with no vegetation, so it can be concluded that vegetation controls the slope stability in the area as also confirmed in previous studies 8,12,28,31 . Updated landslide inventory. In this study, most of the existing landslides which were mapped previously have been also detected in InSAR analysis. Some new landslides possessing high deformation velocity (V slope ) were mapped and several landslides' boundaries were modified according to PS points and validated during fieldwork. The updated landslide inventory was classified into two classes according to deformation rate. Although most landslides have displacement to some extent (Figs. 13, 14), some zones are deforming at the rate of > 25 mm/yr. Those landslides which have high deformation velocity were included in the potential class and other landslides have low deformation values and have information previously classified into existing landslides (Figs. 13, 14). All these potential deforming slopes are complex having different kinds of downfall movement of material and these sites need further analysis by geoscientists to cope with the Ataabad-like disaster in the future.
Landslide susceptibility mapping. The updated landslide inventory map and selected conditioning factors were incorporated with the frequency ratio algorithm to generate the susceptibility model for the area (Fig. 15). The accuracy of the model was 86.6 percent calculated by the AUC method (Fig. 15). The susceptibility map portrayed that, Hunza-Nagar and Gojal-Passu areas are the most hazardous for landslide activities. The susceptibility map was classified into five classes from very high to low susceptible zones (Fig. 15).

Discussion
In this research work, an examination of a Sentinel-1 C-band data set for the potential contribution of PS deformation maps to landslides along the KKH was assessed from Jan. 2018 to Jan 2022. The deformation velocity along the slope (V slope ) map was incorporated to map a number of new potential landslides and adjustment of limits on documented landslides (Fig. 5). Optical remote sensing data interpretation, existing landslide inventory data, and a fieldwork survey of the area were employed for validation, classification, and further analysis of inventory 12,13,28 . The previous studies done by applying InSAR were the SBAS method in the area by 12 using a 1-year temporal period, and 3 also incorporated the same technique for landslide susceptibility optimization in the China section of KKH. In this study, the PSInSAR approach was tested over 4 years of a temporal period which showed promising results, and the relationship between the mapped landslide body and PSInSAR data reveals a good correlation in the area. In this study, several new landslides were mapped that were missed in previous studies like the Budalas landslide having a mean > 50 mm/yr. These updates offer essential information to end-users and stakeholders for the proper planning of risk mitigation measures because greater PS densities are returned in areas with an urbanized environment and a well-developed road network.
Multi-temporal interferometry has many advantages to studying landslides but also encounter some limitation in different nature of areas. Landslides frequently happen in difficult environmental conditions for temporal InSAR applications (e.g. vegetated slopes, steep and rough topography), where InSAR analysis faces significant uncertainty in estimates of ground motions 44 . This limitation can be solved by the integration of long-wavelength ALOS/PALSAR-2 SAR data for temporal decorrelation in vegetation-covered areas 45 . In this study, it was also assessed that the interferometric processing method is also disabled to detect fast-moving landslides like rockfalls. www.nature.com/scientificreports/ So, we propose that detailed geohazard monitoring and identification along KKH be accomplished using a variety of approaches and extensive datasets. Conditioning factors analysis is significant to study landslides. In this research work, the six most responsible factors were evaluated to figure out the correlation with slope failure activities in the area according to all previous studies 7,8,12,32 . In the low terrains, debris flow events are found in many locations along KKH 46 , and 15-60 degrees of sloppy traces were traced to have a more destabilize slope 32 . During rainy seasons loose material from historical landslides flows down and damages KKH frequently. The unstable sections of old landslides in gentle to steep slopes move down slowly found in InSAR results having a huge number of PS points on the historically documented landslides. The barren class was assessed high susceptible to landslides 8 , and more than 60% of the terrains along KKH are barren which weathered with time under the sun, rainfall, and other processes and are prone to landslides. Active tectonic nature accounts for several deformation activities along the fault lines 9 . reported several landslides along the MMT. Several potential deformation events were mapped in the vicinity of the fault lines in this study (Figs. 13,14), which show faults as a main trigger. KKH is passing through most seismic active part of the world where fault and shear zones show a strong influence on landslide activities in the area 7 .Quaternary Deposits (Q) and Southern Karakoram metamorphic Complex (SKm) lithological units were estimated as more prone to landslide in this work which was also assessed by 8,32 . The area is semiarid but in the monsoon season high precipitation is experienced by the region and debris flow, scree slope failure and sometime rockfall events damage the road every year.
Although there are multiple advantages to using RADAR remote sensing for landslide detection and mapping, but some limitations also exist. Data must be evenly spaced in time (regular sampling). Although InSAR provides a constant revisiting period, it is not always possible for the InSAR time series to provide a constant time interval since it is normal for the loss or omission of some images from processing, which affect the results. The wavelength of the SAR sensor being utilized limits the rate of deformation that can be observed. Catastrophic landslides that are moving quickly may deform at rates considerably above meters per minute. The current generation of SAR satellites cannot directly monitor this extremely fast rate of motion 47 . The abundance of vegetation is another drawback. More vegetative cover causes volumetric decorrelation to rise, which reduces coherence 48 . Coherence could be raised by adding corner reflectors to the steepest or most dangerous vegetated slopes, employing longer wavelengths like L-band that can pass through vegetation, and/or raising the temporal resolution but it also encounters atmospheric delay problems. This could result in inaccurate measures of deformation without the application of sophisticated atmospheric mitigation methods. For this study area, susceptibility maps were forwarded by some researchers previously purely based on quantitative methods, that are dependent on visually interpreted inventory data. In this kind of approach, maximum chances of missing inventory data are possible. To overcome this limitation and to map all kinds of landslides for the development of a complete landslide inventory, optical remote sensing interpretation techniques, fieldwork, and InSAR application for landslide detection were applied and generated a highly accurate landslide susceptibility map for the area.

Conclusion
This study presents the development of updated landslide inventory and deformation velocity (V slope ) estimation in a 10 km buffer along the KKH from Gilgit to the Khunjerab section in Northern Pakistan. The results of this investigation show that persistent Scatterer Interferometry (PSInSAR) can significantly update landslide inventories. This InSAR advance works efficiently to modify landslide boundaries, assess their state of activity, and to better understand sliding kinematics. The deformation velocity (V slope ) was absorbed in 364 mm/yr highest in the area. The Budalas area, Khanabad region, Mayoon landslide, and Attabad area are the highly deforming zones that need to be investigated and mitigate future disasters. Landslides possessing displacement > 25 mm/ yr were considered high risk and 29 landslides were mapped and redefined in this study. The PSInSAR-based updated landslide inventory and triggering factor were formulated through the FR model to generate a susceptibility map for the area, which was classified into five zones from very high to low susceptible. Slope, SKm, Q, and Pm lithological units, bareness, and the seismic zones along the fault lines are the most responsible parameters  Figure 12. Graphical representation of Table 2. Those classes which possess zero % landslides were ignored as you can see in the classes' column in Table 2.